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Abstract 

By means of contact dynamics simulations, we investigate a dense packing composed 
of polyhedral particles under quasistatic shearing. The effect of particle shape is an- 
alyzed by comparing the polyhedra packing with a packing of similar characteristics 
except for the spherical shape of the particles. The polyhedra packing shows higher 
shear stress and dilatancy but similar stress-dilatancy relation compared to the 
sphere packing. A harmonic approximation of granular fabric is presented in terms 
of branch vectors (connecting particle centers) and contact force components along 
and perpendicular to the branch vectors. It is found that the origin of enhanced 
shear strength of the polyhedra packing lies in its higher force anisotropy with 
respect to the sphere packing which has a higher fabric anisotropy. Various con- 
tact types (face-vertex, face-face, etc) contribute differently to force transmission 
and fabric anisotropy. In particular, most face-face contacts belong to strong force 
chains along the major principal stress direction whereas vertex-face contacts are 
correlated with weak forces and oriented on average along the minor principal stress 
direction in steady shearing. 
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1 Introduction 



Many recent numerical studies of granular media deal with model systems 
composed of spherical particles. The use of simplified particle shapes and 
contact interactions is needed in order to focus on the collective behavior 
of particles which is at the origin of many specific properties of granular 
materials. On the other hand, the numerical treatment of complex particle 
shapes by discrete element methods was until very recently out of reach due 
to demanding computational resources. There is presently, however, consider- 
able scope for the numerical investigation of complex granular packings. This 
is not only due to available computer power and memory but also because 
during more than two decades of intense research in this field, many funda- 
mental aspects of granular media have already been established for simpli- 
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been analyzed for circular particles (in 2D) and spheres (in 3D). Hence, a recur- 
rent issue today is how robust these findings are with respect to particle shape 
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The issue of shape effect opens actually the door to a vast and substan- 
tial scientific domain given a multitude of potential particle morphologies. 
Several well-known examples are elongated and platy shapes (occurring in 
biomaterials and pharmaceutical applications), angular and facetted shapes 
(occurring in geomaterials) and nonconvex shapes (occurring in sintered pow- 
ders). The macroscopic shear behavior is considerably influenced by parti- 
cle shape. Rounded particles enhance flowability whereas angular shape is 
susceptible to improve s hear streneth, a factor of vit al importance to civil- 
engineering applications (INouguier-Lehon et al.l 20031 ] ). In many engineering 
applications the particle shapes need to be optimized i n order to increase 
jerformance (Marklandl |l981 |. Wu and Thompson 2000l | jLim and MacDowel 
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In this paper, we employ the contact dynamics method to investigate the 
slow shear behavior of granular media composed of polyhedral particles. The 
facetted shapes give rise to a rich microstructure where the particles touch 
at their faces, edges and vertices. Face-face contacts are expected to play a 
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major role in force transmission and statics of polyhedra by accommodating 
long force chains that are basically unstable in a packing composed of spheres. 
In order to isolate the effects arising from particle shape, the data from the 
polyhedra packing will be compared with a packing of spherical particles that, 
apart from the particle shape, is identical in all respects (preparation, friction 
coefficients, particle size distribution) to the polyhedra packing. Both packings 
are subjected to monotonous triaxial compression. 

The numerical procedures will be presented with a brief technical introduction 
to the detection and treatment of contacts between polyhedra in the frame- 
work of the contact dynamics method. We will consider the stress-strain and 
volume-change behavior. The harmonic approximation of the fabric and addi- 
tive decomposition of the stress tensor into fabric and force anisotropies will 
be presented in detail. This will allow us to assess in clear terms the origins of 
shear strength in the polyhedra packing from fabric and force anisotropies in 
comparison to the sphere packing. The probability density functions of normal 
forces will be studied and compared between the two assemblies. Finally, we 
will focus on the contact networks of polyhedral particles and the role played 
by different contact categories with respect to force transmission. 



2 Numerical method 



In this section we briefly introduce the contact dynamics (CD) method with 
polyhedral particles and the numerical procedures used for sample prepara- 
tion. 



2. 1 Contact dynamics method with polyhedra 



The CD method is based on implicit time integration and no nsmooth formula- 



tion o f mutual exclus i on and dry friction betw e en particles (IJean and Moreau 



1992] . lMoreaul [1994j . iRadiai and Rouxl [1999j |. iDubois and Jeanl [2003J ). The 



equations of motion are form ulated as differ ential inclusions in which velocity 
jumps replace accelerations ( Moreau 19941 ] ) . The unilateral contact interac- 
tions and Coulomb friction law are represented as set-valued force laws. The 
implementation of the time-stepping scheme requires the geometrical descrip- 
tion of each potential contact in terms of contact position and its normal unit 
vector. 



At each time step, all kinematic constraints implied by enduring contacts are 
simultaneously taken into account together with the equations of motion in 
order to determine all velocities and contact forces in the system. This problem 
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is solved by an iterative process pertaining to the non-linear Gauss-Seidel 
method that consists of solving a single contact problem, with other contact 
forces being treated as known, and iteratively updating the forces until a given 
convergence criterion is achieved. The method is thus able to deal properly 
with the nonlocal character of the momentum transfers resulting from the 
impenetrability of the particles and friction law. 



The CD method is unconditionally stable due to its inherent implicit time 
integration method. The uniqueness of the solution at each time step is not 
guaranteed for perfectly rigid particles. However, by initializing each step with 
the forces calculated in the preceding step, the variability of admissible so- 
lutions shrinks to the numerical resolution. In the discrete element methods 
based on molecular dynamics, this "force history" is, by construction, included 
in the particle positions. 



The treatment of a contact interaction between two particles requires the 
identification of the contact zone and a "common plane". For rigid parti- 
cles it is possible to define this contact zone by a finite set of points. Before 
applying the contact detection algorithm between a pair of particles of irreg- 
ular shapes, a "bounding box" method is used to compute a list of particle 
pairs potentially in contact. Then, for each pair, the first step is to determine 
i f an over l ap ex i sts through a 3D extens ion of the "shadow overlap method" 
( iSaussind [2004| , iDubois and Jeanl [2003| ) . Several algorithms exist for overlap 



determination between convex pol v hedra (jCundall and Strackl | 1979| .ICundall 
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baussme et al.1 [20061 ]. iPeraled [2007]). When an overlap occurs, the contact 



plane is determined by computing the intersection between the two particles. 



The contacts between polyhedral particles belong to different categories, namely 
face-face, edge-face, vertex-face, edge-edge, vertex-vertex, vertex-edge; see Fig. 
[TJ The vertex- vertex and vertex-edge contacts are practically absent. In all 
cases, we determine one, two or three contact points which provide a good 
description of the contact zone. In this paper, the vertex-edge and edge-edge 
contacts are referred to as "simple" contacts whereas the edge-face and face- 
face contacts are treated as "double" and "triple" contacts since their repre- 
sentation involves 2 and 3 distinct points on the common plane, respectively. 



For our simulations, we used the LMGC90 which is a multipurpose software 
developed in Montpellier, capable of modeling a collection of deformable or 
undeformable par t icles of various shapes by different resolution algorithms 
( IDubois and Jean] [2003J). 
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Fig. 1. Different types of contacts between two polyhedra. 
2.2 Sample preparation 

We generate two numerical samples. The first sample (SI) is composed of 
36933 polyhedra. The particle shape are taken from a library of 1000 digitalised 
ballast grains provided by the French Railway Company SNCF. Each particle 
has at most 70 faces and 37 vertices and at least 12 faces and 8 vertices. Fig. 
[2] shows several examples of the polyhedral particles used in the simulations. 
The size of a particle is defined as two times the largest distance between the 
bary center and the vertices of the particle, to which we will refer as "diameter" 
below. We used the following size distribution: 50% of diameter d m i n = 2.5 
cm, 34% of diameter 3.75 cm, 16% of diameter d max = 5 cm. This distribution 
represents an approximation of that of railway ballast grains. The sample 
contains 7.1 10 5 vertices and more than 10 6 faces, the average numbers being 
20 and 35, respectively. The second sample (S2) is composed of 19998 spheres 
with exactly the same size distribution as in SI. Fig. [3] shows a snapshot of the 
two samples in equilibrium state after deposition and isotropic compression 
under a constant stress of <To = 10 4 Pa in a rectangular box at zero gravity. 

The coefficient of friction is 0.5 between the particles in both samples and 
with the walls. The normal and tangential coefficients of restitution are 
0. The zero restitution simplifies the deposition and compaction process by 
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Fig. 2. Examples of polyhedra used in the simulations. 




Fig. 3. Snapshots of the two packings SI (polyhedra) and S2 (spheres). The walls 
are not shown 

enhancing dissipation during dynamics rearrangements. The initial value of 
the solid fraction is p ~ 0.6 in both samples. Both samples have a nearly 
square bottom of side such that L « I and an aspect ratio H/L ~ 2, where H 
is the height. The initial configuration is denned by Hq ~ 30-Da/ for SI and 
for S2 with Dm the mean diameter. 

The isotropic samples are subjected to vertical compression by imposing a 
constant downward velocity of 10 cm/s on the upper wall and a constant 
confining stress Oi = a 3 = cr on the lateral walls. Each simulation is stopped 
for a vertical deformation of 30%. The time step was 2.10 -4 s. The CPU 
time was 2 10~ 3 s for SI and 1 10~ 3 s for S2, per particle and per time 
step on an Apple G5 computer. The deformation process can be considered 
to be quasistatic in view of the weak kinetic energy injected into the samples 
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compared to the static pressure. This can be expre ssed more generally through 
the inertial number defined as (jGDR-MiDil |2004j |): 




where e = H/H is the vertical strain rate, m is the total mass, p is the mean 
pressure and d is the mean particle diameter. In our simulations, we have 
I ~ 10~ 3 , corresponding to the quasistatic limit. 



3 Stress-strain behavior 



In this section, we compare the stress-strain and volume-change behavior be- 
tween the packings of polyhedra (packing SI) and spheres (packing S2). The 
stress and strain variables are defined from numerical data. For the estima- 
tion of the stress tensor, we use the 'Sensorial mome nt" M % of each particle i 
defined by (IMoreaul [19971 ] , IStaron and Radjai 



where f£ is the a component of the force exerted on particle i at the contact 
is the (3 component of the position vector of the same contact c, and 



c, r 



13 

the summation runs over all contact neighbors of particle i (noted briefly by 
c e i). 



It can be shown that the tensorial moment of a collectio n of rigid part icles is 
the sum of the tensorial moments of individual particles (IMoreaul 199710 . Th e 
stress tensor a for a pack ing of volume V is simply given by (IMoreaul 1997 
Staron and Radjail [20051 ]): 



(3) 



where £ c is the branch vector joining the centers of the two touching particles 
at the contact c. Remark that the first summation runs over all particles 
whereas the second summation involves the contacts, each contact appearing 
only once. 

Under triaxial conditions with vertical compression, we have a% > a 2 = 03, 
where the a a are the stress principal values. Using the Cambridge representa- 
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tion, we define the mean stress p and stress deviator q by (lAirey and Wood 
1983) : 



P= + 0"2 + (T 3 ), (4) 

Q= 3(0-1 -0-3)- (5) 

For our system of perfectly rigid particles, the stress state is characterized by 
the mean stress p and the normalized shear stress q/p. 

The cumulative strain components e a are defined by 



f dH' f AH\ 



Wo 
L 



if' \ H, 



10 



. dV , / AL\ 



io 



-0 



e 3 = / _=ln 1 + — , (8) 







where ifo^ an d L are the initial height, width and length of the simulation 
box, respectively and AH = H — H , Al = l — I and AL = L — L are the 
corresponding cumulative displacements. The volumetric strain is given by 

V fdV / AV\ 
e P = I — = !*(! + --), (9) 



where Vo is the initial volume and AV = V — Vq is the total volume change. 
The cumulative shear strain is defined by 

e q = e l -e 2 . (10) 



Figure H] displays the evolution of q/p for the packings SI and S2 as a func- 
tion of e q . For both packings, we observe a classical behavior characterized by 
a hardening behavior followed by (slight) softenin g and a stress plat e au co r- 
responding to the critical state of soil mechanics ([Mitchell and Sogal 20051 ] ). 
The critical-state strength in the case of polyhedra (~ 0.46) is twice as high 
as that of spheres (~ 0.23). This implies that the polyhedra packing has a 
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Fig. 4. The strength parameter q/p as a function of shear strain e q for the polyhedra 
packing SI and sphere packing S2. 
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Fie. 5. The volume change function of shear strain e q for the packings SI et 
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At the critical state, we have ip = <fo = 34° for SI and fo = 18° for S2. 

Figure O shows the volumetric strain e p as a function of shear strain e q in SI 
and S2. In both packings, we observe an early compaction slightly larger in S2 
than in SI. The subsequent dilation is lower in S2 and the critical state with 
isochoric deformation is reached at e q = 0.3. Dilation in SI continues with a 
decreasing rate of volume change but the isochoric plateau is not fully reached. 
The dilatancy can be expressed in terms of the dilation angle if defined by 



sin^ 



:i2) 



We have ip ~ 5° for SI and ip ~ 2.5° for S2 at the stress peak state. 
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Fig. 6. The stress-dilatancy diagram representing the relation between the internal 
angle of friction and the dilation angle for polyhedra and spheres. 



The variation of ip versus (p, a sort of stress-dilatancy diagram (1 Wood! 19901 ] ). 
is displayed in Fig. O for polyhedra and spheres. For both packings, we have 



(13) 



where k is a constant slightly smaller than 1 in both packings. This corre- 
lation between dilatancy and shear stress during stress-strain transients is a 
consequence of energy balance. The mechanical work performed on the sys- 
tem is parti ally dissipated in contac t interactions and partially used in vol- 
ume change (IRadjai and Rouxl 20041 ] ). Several stress-dilatancy relation s have 
been p roposed as flow rules for plastic deformations of granular media (IWood 
1990]). The relation (|T3|) associates the peak state to the largest positive value 
of dilatancy and the critical state to zero dilatancy. It shows the "non asso- 
ciated" character of the flow rule in granular media (an associated flow rule 
implying (p — ijj). 



4 Harmonic representation of the fabric 



The expression of stress tensor in Eq. ([3]) is an arithmetic mean involving 
the branch vectors and contact forces. Hence, in order to analyze the shear 
strength properties of the polyhedra packing compared to the sphere packing, 
we need a statistical description of the granular microstructure (texture or 
fabric) and force transmission. 



In the presence of steric excl usions, the granular microstructure i s highly dis- 
ordered at the particle scale ( Troadec 2002], Troadec et al. 20021 ] ). Since me- 
chanical interactions are governed by contact and friction, the relevant de- 
scriptors of the microstructure are related to the contact network. At the 
lowest order, the contact network is characterized by the coordination num- 
ber z which describes the compactness of a packing. This is a crude scalar 
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Fig. 7. Evolution of the coordination number z as a function of the cumulative shear 
strain e q for polyhedra (SI) and spheres (S2) 
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Fig. 8. The connectivity P{c) of the contact network for the packings SI and S2. 

information in view of the complex arrangement of the particles, but it is 
well-known that the compactness, generally expressed in terms of the solid 
fraction, controls the stress-strain behavior under monotonous shearing. Let 
us remark here that double and triple contact types (see section [2]) are counted 
as single contacts for the coordination number although they are represented 
by two and three contact points, respectively, in the numerical treatment of 
interactions between polyhedra. 

The evolution of z for polyhedra and spheres is shown in Fig. [7] as a function 
of e q . It is remarkable that z is nearly constant in spite of the overall dilation 
in both packings. We have z ~ 5.5 for polyhedra and z ~ 4 for spheres. The 
connectivity of the contact network can be characterized in more detail by the 
fraction P(c) of particles with exactly c contact neighbors. The coordination 
number is the mean value of c : z = J2c c P( c )- Fig. [8] shows P(c) for SI and S2 
in the critical state. The distribution is broader in SI than in S2. This shows 
the wider range of potential equilibrium states in the polyhedra packing. For 
both packings, we observe a peak centered on c = 4 with a higher probability 
forS2. 



Since the shear stress corresponds to the deviation of stress components from 
the mean stress p along different space directions, the coordination number 
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Fig. 9. Geometry of a contact between two polyhedra. 

z as a scalar quantity cannot account for the shear stress and its evolution 
with strain. Indeed, the expression of the stress tensor suggests that the useful 
information for the analysis of shear stress is the density and average force as 
a function of cont act orientation. Such functions c an be expanded in spherical 



harmonics in 3D (jOuadfel and Rothenburgi 20011 ]) 



Let n be the unit vector along the branch vector £ ; Fig. [91 We set 

i = in, (14) 



where I is the length of the branch vector. We remark that the unit vector n 
does not coincide with the contact normal except in the case of spheres. We 
consider the components of the contact force in a local frame defined by n 
and an orthoradial unit vector t: 

f = f n n + f t t, (15) 



where f n and f t are the radial and orthoradial components of the contact 
force, respectively. The writing of Eq. ffTBT) assumes that t is oriented along 
the orthoradial force. 

We now define the angular averages associated with the branch vectors £ and 
contact force vectors /. Let A(Q) be the set of branch vectors pointing in the 
direction Q = (6, <ft) up to a solid angle dQ and N C (Q) its cardinal. The angles 
9 and </> are shown in Fig. [TUJ The angular averages are defined as follows: 



Pn{£l) = ^^ } (16) 



• c 
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Fig. 10. Spherical coordinates. 
(/»>(")= AFToT E (18) 

(/*>(«) = 1v7oT E (I 9 ) 

where N c = J N c (Q)dfl is the total number of contacts, and £ c , /°, and / t c are 
the actual values of branch vector length, radial force and orthoradial force 
for contact c, respectively. 



Under the axisymmetric conditions of our simulations, the four functions de- 
fined in Eq. (fT9l) are independent of </>. Fig. [TT] displays a polar representation of 
these functions in the #-plane for polyhedra (SI) and spheres (S2) at e q = 0.3. 
We observe an anisotropic behavior in all cases except in (£)(9) for S2. A weak 
anisotropy can be seen for SI in the latter case. The peak values occur along 
the compression axis except for (ft) (9) in which the peaks are inclined at 7r/4 
with respect to the vertical. The magnitude of anisotropy is larger for poly- 
hedra compared to spheres except for Pq{9) which is weakly anisotropic for 
polyhedra. 

The simple shapes of the above functions suggest that harmonic approximation 
based on spherical harmonics at leading terms captures their anisotropies. 
There are 9 second-order basis functions Yj^(9, <fi). But only the functions 
compatible with the symmetries of the problem, namely independent with 
respect to and 7r-periodic as a function of 9, are admissible. For Pq{9) 
as a scalar, and (t){9) and (f„){9) as radial components of the vectors, the 
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Fig. 11. Polar representation of density probability function Pq(6), (f n )(9), (ft)(9) 
and (i)(0) for SI et S2 in residual state. 

only admissible functions are Y" ° = 1 and Y" 2 = 3cos 2 # — 1. For (ft) {9) as 
orthoradial component of a vector, the only function independent of <fi and 
perpendicular to Y" ° = 1 and Y" 2 = 3cos 2 # — 1 is sin 29. Hence, within the 
harmonic model of fabric and force, we have 



Pn{6) = —{ l + a [3 cos 2 9- 1] }, 

47T 

(£)(6)=l { l + a t [3 cos 2 9- 1] } 
(fn)(0) = f o {l + a n [3cos 2 £-l] }, 
(ft)(9) = f a t sin2[0], 



(20) 

(21) 
(22) 
(23) 



where and a t are the anisotropy parameters, lo is the mean branch 

vector length, and /o the mean force. The probability density function Pn(6) 
is normalized to 1 (f s Pn(Q)dQ = 1, where S is a sphere of unit radius). The 
values of the anisotropies a, a;, a n and a t can be calculated from generalized 
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Fig. 12. Evolution of anisotropies a, ai, a n and at with e q for packings SI and S2. 
fabric tensors introduced in Appendix IA1 

The evolution of the anisotropies with displayed in Fig. [12] for our 

packings SI and S2. The fabric orientation anisotropy a increases with e q 
and relaxes to a plateau after passing by a pronounced peak. Its value is 
systematically larger for spheres than for polyhedra (by a factor 3 in the 
critical state). The branch vector anisotropy a/ is quite low compared to other 
anisotropies and its value all along shearing is negligible for spheres. It is 
remarkable that a/ for polyhedra declines ( Fig [5]) at the beginning 

of shearing. The radial force anisotropy a n increases as the fabric anisotropy 
and tends to a plateau. But, in contrast to fabric anisotropy, its value is 
higher for polyhedra than for spheres. In other words, the aptitude of the 
polyhedra packing to develop large force anisotropy is correlated with particle 
shape rather than with fabric anisotropy (see section [7]). The orthoradial force 
anisotropy a t has a similar behavior except that it takes considerably higher 
values in the case of polyhedra compared to spheres. In the following section, 
we study the relationship between the fabric and force anisotropies. 



5 Origins of shear stress 



In this section, we analyze the stresses in the framework of the harmonic 
approximation of granular microstructure introduced in the last section. Since 
this representation involves continuous functions of contact orientations, we 
need to express the stress tensor in integral form. The stress tensor as defined 
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in Eq. is an average: 
o af 3 = n c {£ k J$) k , 



(24) 



where n c = N c /V is the number density of contacts, £^ is the a component of 
the branch vector at contact k and £^ is the /3 component of the force vector 
at contact k. The average is taken over all contacts k in the control volume V. 
To express this mean as an integral, we intro duce the joint probability density 



Ppjf.in, f , 1) of the force and branch vectors (IBathurst and Rothenbura [1988 



Rothenburg and Bathurst 1989 ]. Ouadfel and Rothenburgi 2001 1). Then, from 
Eq. (J24J), we have 



v a p = n c J Pn fi (n, f,£) i(n) fp{n,£) n a dtt dfdt, 



(25) 



where dtt = sin 



Equation (125]) can be simplified by integrating out the contribution of I. As- 
suming that / is independent of I (an assumption which is verified with a 
good approximation), we get 



a a/3 = n c J Pnf{n,f) (£){n) fp(n) n a dVt df, 



(26) 



where (£) P nf = J P 



Finally, integration of f[26|) over force vector yields the following expression for 
the stress tensor: 



cr a/ 3 = n c P n (n) (£){n) n a (/g)(n) dtt, 



(27) 



where (f)Pn — I Pnf df. By introducing the average force components (/„) 
and (f t ) in this equation, we get 



(J a p = n c / P n (n) (£)(n) { (f n ){n) n p + (f t )(n) t p }dtt. 



(2f 



This writing of the stress tensor involves the functions previously introduced 
with the harmonic representation of the fabric (Eqs. (|21|) . f[2"2"j) . (|2"51) and f|2^]) ). 
Inserting these functions in the integral expression (Eq. [281) and given the 
definitions of mean stress p and stress deviator q in Eq. (J5j), one gets 



V - nJofo 
q 2 



P 



a + a/ + a n + at), 



(29) 
(30) 
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Fig. 13. The normalized shear stress q/p as a function of shear strain e q for the 
packings SI and S2 both from direct simulation data and theoretical prediction of 
Eq. ([3D]). 



where the cross products (aai, aa n and aa t ) among the anisotropies have 
been neglected. Our simulation data are in quantitative agreement with this 
"stress-force-fabric" relation (a t e rm coined by Rothenburg and B athurst in 
( Bathurst and Rothenburg 1988 ]. Rothenburg and Bathurst 1989j |) ) both for 
spheres and polyhedra, as shown in Fig. Q21 all along the shear. We note that 
the theoretical fit would have been less satisfactory for polyhedra if the branch 
vector length anisotropy a; were omitted from the description. 



Equation fl30l is interesting as it exhibits the two origins of shear stress in 
a granular system: 1) the fabric anisotropies a and ai, related to the branch 
vector and 2) the force anisotropies a n and a t , related to the contact force. 
Figure d2] shows that the values of these anisotropy parameters underlying 
the shear stress depend on the particle shape. In particular, the total force 
anisotropy a n + a t compared to the total fabric anisotropy a + a; is much 
higher in the case of polyhedra. In the critical state, we have a n + at — 0.88 
and a+ai ~ 0.2 for polyhedra, a n + a t — 0.26 and a+ai ~ 0.24 for spheres. The 
high value of the force anisotropy in the case of polyhedra comes from both 
radial and orthoradial components whereas in the sphere packing at — 0.05 
is much less important than a n ~ 0.21. This suggests that friction is more 
directly involved in force transmission in the polyhedral packing than in the 
sphere packing. The strong contribution of force anisotropy to the polyhedra 
packing is a particle shape effect related to the face-face contacts which carry 
most strong forces. This point will be analyzed in more detail below. 



6 Force distributions 



In this section, we study the probability density functions (pdf's) P(f n ) f° r 
sphere and polyhedra packings. Fig. [TH shows typical maps of normal forces 
in a portion of both packings in the critical state. The 3D force chains can 
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Fig. 14. Force maps in a portion of the packings SI (right) and S2 (left). The 
segments are branch vectors with thickness proportional to the normal force, and 
gray level proportional to the depth of field. 
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Fig. 15. Probability density functions of normal forces in the packings of spheres 
and polyhedra. 

be observed in both packings, but they seem more tortuous in the case of 
polyhedra. 

The normal force pdf 's are shown in Fig. [15] on log-linear and log- log scales at 
e q = 0, 3. In both pdf 's, the strong forces, i.e. forces above the mean normal 
force (/„), fall off exponentially: P(f n ) oc e~^ n f^ n \ with (3 ~ 0.9 for SI 
and j3 ~ 1.1 for S2. In contrast, the shapes of the pdf's in the range of 
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weak forces (/„ < (f n )) are radically different. In the sphere packing, the pdf 
slightly bends down as f n — > but does not tend to zero. We observe also 
a small peak close to the mean force. This is consistent with s everal other 
nume ri cal and experimental obs e rvations for isotropic packings (ILoyol et al. 



1999 



2004 



Bardenhagen et al. 2000 



Majmudar and Behringer 



Antony! [200 lj . ISilbert et alJ |2002| . iMetzgei 
2005) ) . In the case of polyhedra, the number 



of weak forces bends up as the force tends to zero. For both packings, the range 
of weak forces is well approximated by a power-law distribution : 



P(fn) OC 



(31) 



with a = —0.24 for SI and a = 0.05 for S2. The divergence of the number of 
weak forces in SI should be attributed to the polyhedral shape of the particles 
favoring the arching effect an hence a higher fraction of weak forces. The 
coefficient of friction has a similar effect though to a lesser extent. We find, 
however, that in both systems the fraction of weak forces (/„ < (f n )) is about 
60%. 



7 Contact networks of polyhedral particles 



In the case of the polyhedra packing, it is interesting to investigate the orga- 
nization of the contact network in terms of simple, double and triple contacts. 
The respective fractions of these contact types and their contributions to the 
structural anisotropy and force transmission are the key quantities for un- 
derstanding the effect of particle shape on the shear strength properties of 
granular media. In fact, one expects that the triple (face-to-face) contacts 
play an essential role in force transmission. This f eature was observed in the 
case of polygon packings for side-to-side contacts f Azema et al. [20071 ] ). 



Considering the discrete expression of the stress tensor in Eq. (jHJ) and re- 
stricting the summation to each contact type allows us to perform an additive 
decomposition: 

<r = cr s + cr d + a,. (32) 



where the subscripts s, d and t design the respective contributions of simple, 
double and triple contacts. The corresponding stress deviators q s , qd and q t 
are then calculated and normalized by the mean stress p. Fig. [TH] shows the 
evolution of partial shear stresses q s /p, qd/p and q t /p as a function of shear 
strain e q . The contribution of simple contacts is larger for double and triple 
contacts. However, the double and triple contacts support together the largest 
portion of the overall shear stress, i.e. qd + qt > q s - 
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Fig. 16. Evolution of partial shear stresses as a function of shear strain for simple 
(s), double (d) and triple (t) contacts, as well as the total shear stress (s+d+t). 
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Fig. 17. Proportions k s , kd and k t of simple, double and triple contacts (dashed 
lines), and the relative average forces f s , fd and f t (full lines) supported by each 
contact type function of shear strain e„. 

The partial shear stress supported by each contact type depends on both the 
number of its contacts and their mean force. Fig. [171 shows the proportions k s , 
ka and k t of simple, double and triple contacts as a function of shear strain. 
k s declines during shear from 0.75 to 0.71 whereas kd and kt increase from 
0.14 to 0.15 and from 0.11 to 0.14, respectively. Hence, the critical state is 
characterized by k s ~ 0.7 et k t ~ kd — 0.15. Fig. [TT1 also shows the relative 
mean forces / s , fd and f t defined by 



f s = k s (f n ) s /(f n ), (33) 
f d = k d (f n ) d /(fn), (34) 
ft = k t (f n )t/(f n ), (35) 

where (f n ) s , (f n )d and {f n )t correspond to the mean normal forces of simple, 
double and triple contacts. We see that f s declines slightly with strain but 
is nearly two times larger than f t and 2.3 times larger than fd in the critical 
state. We have f s ~f d + f t . Hence, the lower contribution of triple contacts 
with respect to shear stress can be attributed to both the low level of the mean 
force (ft < 0.3) sustained by this class and to their weak number (< 15%). 
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Fig. 18. Evolution of the signed anisotropies a' of simple (s), double (d) and triple 
(t) contacts as a function of shear strain e q . 

Following the same procedure as for the stress tensor, we now perform a similar 
decomposition of the fabric tensor F, defined by Eq. (1A.4j) . into three terms: 



where F s , F d and F t are the contributions of simple, double and triple contacts. 
The corresponding anisotropies a s , a d and a t can be extracted, but since the 
principal directions of these partial fabric tensors are not necessarily identical 
to those of the overall fabric tensor, we define the "signed" anisotropies by 
multiplying each partial anisotropy by a phase factor cos 2(6 f — 



Figure [18] shows the evolution of signed anisotropies of the three contact 
classes. We see that a' d and a' t increase with shear strain and tend to the 
limit value ~ 0.04. As to a' s , we observe an initial increase followed by rapid 
decrease and change of sign at e q ~ 0.2. This evolution means that during 
shear the branch vectors of simple contacts tend to become increasingly per- 
pendicular to the major principal direction (the direction of compression). A 
map of contact forces projected along the branch vectors is displayed in Fig. 
[T9l in different colors according to the type of contact. The triple contacts, 
despite their lower proportion, appear clearly to be correlated in the form of 
long chains across the packing. These are mostly parallel to the direction of 
compression. We also observe a large number of weak forces mainly at simple 
contacts. 

The pdf 's of normal forces are shown in Fig. [20] separately for simple, double 
and triple contacts. The three contact types are involved in strong and weak 
networks. The strong forces have in all cases an exponential behavior but a 
major difference is observed in the range of weak forces where the proportion 



F = F S + F d + F u 



(36) 
(37) 
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Fig. 19. Map of contact forces projected along branch vectors at e q = 0.4. Line 
thickness is proportional to the force. The simple, double and triple contacts are in 
red (dark gray), in green (light gray) and in blue (black). 
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Fig. 20. Probability distribution functions of radial forces at simple (s), double (d) 
and triple (t) contacts on log-linear (a) and log-log (b) scales. 



of simple contacts prevails. This correlation between simple and weak contacts 
is interesting as it clearly reveals the contrast between simple contacts, on one 
hand, and double and triple contacts, on the other hand, in the organization 
of the force network. 
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Fig. 21. Proportions kg , kf and /cf of simple (s), double (d) and triple (t) contacts 
in the strong network (S) and the corresponding proportions kY , kY and kY in 
the weak network (W) function of shear strain. 

In order to situate the simple, double and triple contacts with respect to the 
force network, we have plotted in Fig. [21] the proportions kg, kf et kf of the 
three contact sets in the strong network and the corresponding proportions 
kY , kY et kY in the weak network as a function of shear strain e q . It is 
interesting to note that the proportion of weak simple contacts is quite high 
(~ 0.55). The proportions kY et A;f of weak and strong double contacts are 
identical (~ 0.07). Finally, we see that most double contacts belong to the 
strong network (kf ~ 2kf )■ 



8 Conclusion 



In this paper, granular materials composed of irregular polyhedral particles 
were numerically investigated. Macroscopic and microstructural properties 
were analyzed by (1) direct comparison with a similar packing composed 
of spherical particles and (2) characterization of contact networks and force 
transmission. A novel finding of this work is that the origin of enhanced shear 
strength in a polyhedra packing compared to a sphere packing lies in force 
anisotropy induced by particle shape. The fabric anisotropy associated with 
the network of branch vectors is lower in the polyhedra packing. This find- 
ing extends the results of a p revious study of pe ntagonal particles in two 



dimensions to three dimensions Azema et al. 2007 . In other words, the force 



anisotropy, partially underlying shear strength, is mainly controlled by the 
fabric anisotropy in a sphere packing. This mechanism breaks down to some 
extent in a packing of polyhedra where force anisotropy results mainly from 
the "facetted" particle shape. 

The face-face contacts were shown to belong mostly to the strong force net- 
work. The local equilibrium structures involving face-face and edge-face con- 
tacts accommodate force lines that are basically unstable with spherical par- 
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tides. Hence, the term "arching" seems to be more adapted to the description 
of force patterns in an assembly of polyhedra than in an assembly of spheres. 
This effect is crucial for the probability density function of normal forces in 
the range of weak forces that is well approximated by a decreasing power-law 
in the case of polyhedra. 

In this investigation the polyhedra were irregular with a given number of 
faces, edges and vertices. These shape parameters can now be changed and 
the resulting packings can be analyzed along the same lines as in the present 
investigation. Since the face-face contacts seem to play a key role, it would be 
interesting to consider irregular polyhedra with less faces in number but with 
larger areas. From a mechanical point of view, there should be little difference 
between a small face and a vertex. The best shape from the shear strength 
viewpoint can be obtained with a large number of faces of large area, but 
these two conditions can not be realized at the same time. It seems thus that 
an optimal polyhedral shape should exist with a number of faces of not two 
low areas. The work is under way to elucidate this point and other aspects of 
the problem concerned with packing structure by systematically changing the 
particle shape parameters. 

We acknowledge assistance by F. Dubois with the LMGC90 platform em- 
ployed for the simulations, as well as the precious help of V. Richefeu with 3D 
visualization of forces. This work was funded by the French Railway Society, 
the SNCF, and the Region Languedoc-Roussillon of France. 



A Fabric tensors 



The anisotropies a, a n , a t and ai can be calculated from the tensors F, H^ n \ 

H w and defined by (IBathurst and Rothenburgj 1988j |. lRothenburg and Bathurst 
|l989t . lOuadfel and RothenbureJ [200ll ]) : 



F al 3 = J Pq(9) n a n p dtt, 
s 


(A.l) 


H^ = J(f n )(0) n a n p <m, 
s 


(A.2) 


H% = J (ft) (9) n a t p dQ, 
s 


(A.3) 


H% = J{t){B) n a n p dQ. 
s 


(A.4) 
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Using the equations (l2~Tj) . f|22|) . (1231) and (1231 . it is then easy to show that the 
corresponding anisotropies are : 



5 F 3 — F x 

a = — — , A. 5 

2 trF ' v ' 

c rr(") tr(") 

a " = 2 frflW ' (A ' 6) 

a t = d / , (A.7) 

a , = ^ (Ag) 

where trH^ = (/), trF = 1 et tr-FfW = 4 ■ 
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